SU-4240-629 
SCCS-754 
NBI-HE-96-12 
March 1996 



The Flat Phase of Crystalline Membranes 



Mark J. Bowick, Simon M. Catterall 
Marco Falcioni and Gudmar Thorleifsson 

Department of Physics 
Syracuse University 
S yracuse, NY 13244-1130, US A 
http : / / www . phy . syr . edu/| 



Konstantinos N. Anagnostopoulos 

The Niels Bohr Institute 
Blegdamsvej 17 
DK-2100 Copenhagen 0, Denmark 



tittp : //www . nbi . dk/ 



Abstract 

We present the results of a high-statistics Monte Carlo simulation of a phantom 
crystalline (fixed-connectivity) membrane with free boundary. We verify the ex- 
istence of a flat phase by examining lattices of size up to 128 2 . The Hamiltonian 
of the model is the sum of a simple spring pair potential, with no hard-core 
repulsion, and bending energy. The only free parameter is the the bending rigid- 
ity k. In-plane elastic constants are not explicitly introduced. We obtain the 
remarkable result that this simple model dynamically generates the elastic con- 
stants required to stabilise the flat phase. We present measurements of the size 
(Flory) exponent v and the roughness exponent £. We also determine the criti- 
cal exponents r\ and r\ u describing the scale dependence of the bending rigidity 
(K(q) ~ q~ v ) and the induced elastic constants (X(q) ~ fi(q) ~ q Vu ). At bending 
rigidity k = 1.1, we find v = 0.95(5) (HausdorfT dimension da = 2/v = 2.1(1)), 
( = 0.64(2) and r\ u = 0.50(1). These results are consistent with the scaling 
relation ( = (2 + ?? u )/4. The additional scaling relation rj = 2(1 — £) implies 
r] = 0.72(4). A direct measurement of rj from the power-law decay of the 
normal-normal correlation function yields rj w 0.6 on the 128 2 lattice. 
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1 Introduction 



The physics of flexible membranes, two-dimensional surfaces fluctuating in three 
dimensions, is extremely rich, both on the theoretical side, where there is a nice 
interplay between geometry, statistical mechanics and field theory, and on the ex- 
perimental side, where model systems abound. 

The simplest examples of 2-dimensional surfaces are monolayers, or films — these 
are strictly planar systems. The statistical behaviour of monolayers falls into three 
distinct universality classes, depending on the microscopic interactions of the system. 
There are crystalline monolayers, with quasi-long-range translational order and long- 
range orientational order, hexatic monolayers, with short-range translational order 
but quasi-long-range orientational order, and fluid monolayers, with short-range 
translational and orientational order 

Physical membranes, 2-dimensional surfaces fluctuating in a 3-dimensional em- 
bedding space, are expected to have an equally rich phase structure. The simplest 
class of membranes is the crystalline class, in which topological defects, dislocations, 
and disclinations, are forbidden. At the discrete level crystalline membranes may 
be modelled by triangulated surfaces with fixed local connectivity. In this paper we 
will be concerned with the properties of a particularly simple model of a phantom 
(non self-avoiding) crystalline membrane, with emphasis on a critical analysis of the 
existence and stability of an ordered flat phase of the model for sufficiently large 
bending rigidity. 

Crystalline membranes, also known as tethered or polymerised membranes, are 
the natural generalisation of linear polymer chains to two dimensions. Polymer 
chains in a good solvent crumple into a fractal object with Hausdorff dimension 5/3 
(Flory exponent v = 3/5). Crystalline membranes with in-plane solid- like elasticity, 
on the other hand, are predicted to exhibit quite different physical properties from 
their linear counterparts. In particular, they are expected to have a remarkable 
low-temperature ordered phase. This ordered, or flat, phase is characterised by 
long-range order in the orientation of surface normals. At high temperature, or 
equivalently low bending rigidity, phantom crystalline membranes will entropically 
disorder and crumple. Separating these two phases should be a crumpling transition, 
whose precise nature for physical membranes is still not fully understood. 

Inorganic examples of crystalline membranes are thin sheets (< 100 A) of 
graphite oxide (GO) in an aqueous suspension [Q, || and the rag- like structures 
found in MoS 2 i- 

There are also beautiful biological examples of crystalline membranes such as 
the spectrin skeleton of red blood cell membranes. This is a two-dimensional tri- 
angulated network of roug hly 70,000 plaquettesQ. Actin oli gomers form nodes and 
spectrin tetramers form links |y. Crystalline membranes can also be synthetised 
in the laboratory by polymerising amphiphillic monolayers or bilayers. For recent 

lr The 128 2 lattice we simulate has 32,766 plaquettes. 
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reviews see |6|-|8]]. 

The existence of an ordered phase in a two dimensional system is remarkable, 
given the Mermin- Wagner theorem. In fact, if one thinks of surface normals as 
spins, then the membrane bending energy is akin to a Heisenberg interaction, and 
the 2d-Heisenberg model is well known to have no ordered phase. What stabilises 
the flat phase in a crystalline membrane? There are several appealing arguments 
for the existence of a stable ordered phase in crystalline membranes ||. Out-of- 
plane fluctuations of the membrane are coupled to in-plane "phonon" degrees of 
freedom because of the non-vanishing elastic moduli (shear and compressional) of a 
polymerised membrane. Bending of the membrane is inevitably accompanied by an 
internal stretching of the membrane. Integrating out the phonon degrees of freedom 
one finds long-range interactions between Gaussian curvature fluctuations which 
stabilise a flat phase for sufficiently large bare bending rigidity ||. Alternatively 
both the elastic constants and the bending rigidity pick up anomalous dimensions 
in the field theory sense. The bending rigidity receives a stiffening contribution at 
large distances via the phonon coupling. This competes with the usual softening of 
bending rigidity seen in fluid membranes, with the net result being an ultraviolet- 
stable fixed point — the crumpling transition. From the magnetic point of view 
membrane models are constrained spin systems, since the spins must be normal 
to the underlying surface, and the constraints are essential to the stability of the 
ordered phase. 

This viewpoint is supported by self-consistent calculations for the renormali- 
sation of the bending rigidity, by large-d calculations, where d is the dimension of 
the embedding space, and by e-expansion calculations, where e = 4 — D, with the 
internal dimension D of the membrane being 2 for physical polymerised membranes. 

The construction of a discrete formulation of a crystalline membrane is essential 
for numerical simulations, and revealing for comparison with spin systems. In the 
simplest discretised version the membrane is modelled by a regular triangular lattice 
with fixed connectivity, embedded in d-dimensional space. Typically the link-lengths 
of the lattice are allowed some limited fluctuations. This may be modelled by tethers 
between hard spheres or by introducing some confining pair potential with short- 
range repulsion between nodes (monomers). The bending energy is represented by 
a ferromagnetic-like interaction between the normals to nearest-neighbour "plaque- 
ttes" §. 

In the random surface literature, on the other hand, it is more common bto 
consider a simple Gaussian spring model of the pair potential with vanishing equi- 
librium spring length. In this case the minimum link-lengths are unconstrained. A 
priori this model seems to be rather different from those above; one may worry that 
it is pathological in some sense. In particular one may ask if this model can ever 
generate the effective elastic constants which are required to stabilize a flat phase 
at large bending rigidity. To illustrate this concern consider the infinite bending 
rigidity limit of the model. In this limit the system becomes planar — out-of-plane 
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fluctuations are completety suppressed. Since the pair potential allows it, what 
prevents all nodes from collapsing on each other within the plane of the membrane 
itself? The answer can be seen by considering what happens if a node moves across 
a bond, or if two nodes are exchanged. In this case the normals to some triangular 
plaquettes will inevitably be inverted. This generates a prohibitive bending energy 
cost and is effectively forbidden. The model thus dynamically generates a hard-core 
repulsionP], and may be thought of as having an effective equilibrium spring length. 
In Appendix A we show that small fluctuations about a finite microscopic equilib- 
rium spring length a yield elastic constants proportional to (a/e) 2 , where e is the 
intrinsic lattice spacing. For e ~ a these constants are finite. For the model we 
consider, with vanishing microscopic a and finite e, the heuristic arguments above 
suggest one should replace a by an effective equilibrium spring length. 

Monte Carlo simulations have in fact established strong evidence for a continu- 
ous crumpling transition in the model with a = 0. Most of the simulations, however, 
have focused on the crumpling transition itself. Not much effort has been made to 
establish rigorously the existence and the properties of the flat phase. In the rest 
of this paper we present evidence that there is indeed a stable flat phase in this 
remarkably simple model of a crystalline membrane. Furthermore we show that 
the requisite elastic constants are dynamically generated with the correct scaling 
behaviour. 

The paper is organised as follows: in Section ||| we discuss the discrete model, 
the numerical simulations and the choice of the boundary conditions. In Section |3| we 
review the evidence for the crumpling transition. In Section || we discuss the Monge 
representation of a surface and the theoretical predictions for asymptotically flat 
elastic surfaces. We then discuss our numerical results for the roughness exponent, 
the phonon fluctuations and the normal-normal correlation function. In Appendix A 
we give a calculation of the elastic constants of a discrete soft-core Gaussian model. 
In Appendix B we discuss of the methods used to measure the geodesic distance for 
the correlation function and in Appendix C we describe the parallel algorithm used 
to simulate the lattices with largest size. 

2 Model 

To describe a discrete polymerised surface we arrange iV particles (monomers) in 
a regular triangulation of a 2— D manifold (see Fig. |l|a). The 2— D surface is then 
embedded in a d-dimensional space where it is allowed to fluctuate in all directions. 
Each monomer is labelled by a set of intrinsic coordinates a = (cri, 02), with respect 
to a set of orthogonal axes in the 2— D manifold. The position in the d-dimensional 
embedding space is given by the vector x ff . We will treat the case d = 3. 

In general, the Hamiltonian which describes these models has two terms: a pair 

2 We thank one of the referees for informative observations on this point 
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Figure 1: a) The intrinsic connectivity structure of the mesh, b) The labelling 
scheme: <r, p, r, fi are the intrinsic coordinates while a, b label the triangles. 




potential and a bending energy term, 

n = n T + n B . (l) 

A commonly studied pair potential in the literature [|l(| is the tethering potential 
with hard core repulsion, given by 

Ht= £ ydx.-x^l), (2) 

{aa') 

with 

V{R) = l a<R<b ,R=\x a -x a >\ (3) 
I oo i? > b 

where R is the link length, a is the hard-core radius and b is the tethering length. 

We consider, instead, a model where the tethering potential is replaced by a 
simple Gaussian spring potential. The bending energy is the usual ferromagnetic 
interaction between the normals to the faces of the membrane, namely 

H= ^2 \x a - x CT /| 2 + (1 - n a • n 6 ) . (4) 

(era') (ab) 

Here k is the bending rigidity, a and b label the faces of the surface and n a is the unit 
normal to the face a. Both sums extend over nearest neighbours (see Fig. [!]&). This 
action has no explicit short scale cut-off length or hard-core repulsion. Since there 
is no self-avoidance it represents a phantom surface. This action was investigated 
originally by Ambj0rn et al. [11| and is inspired by the Polyakov action for Euclidean 
strings with extrinsic curvature ||l2, 13[. 



The main justification for studying phantom surfaces as models of realistic 
membranes is that in the flat phase self-avoidance is irrelevant JiO|, Although 
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self-avoidance is expected to change the scaling behaviour at the critical point k c 
and the nature of crumpled phase, it does not affect the scaling behaviour in the 
flat phase. As membranes with strong self-avoidance (and the resultant non-local 
interactions) are harder to simulate numerically, it is sensible to leave out self- 
avoidance when possible. 

The partition function of this model is given by the trace of the Boltzmann 
weight over all possible configurations of the embedding variables x: 



Here x cm is the centre of mass, which is held fixed to eliminate the translational 
zero mode. Expectation values are given by 



We will consider the case of surfaces with the topology of the disk and free boundary 
conditions. Most experimental realisations of membranes have either this topology 
or spherical topology. Our choice also has certain technical advantages. To describe 
the flat phase we need to measure the finite-size-scaling of the thickness of the 
surface and the asymptotic behaviour of the normal-normal correlation function. 
With spherical or toroidal topology one would have to subtract the effects of the 
global shape of the surface. 

2.1 Numerical Simulations 

To evaluate the integral of Eq. @ we use the Monte Carlo algorithm with a Metropo- 
lis update. In our case the Metropolis update corresponds to changing the position 
of a node by a trial vector e chosen randomly (and uniformly) in a box of size 
(25) 3 centred on the old node position. The update is accepted if the change in the 
Hamiltonian is such that 



where r is a uniformly distributed random variable with values in the interval [0, 1). 
The value of 5 is adjusted to keep the acceptance ratio around 50%. For a surface 
of linear size L, one Monte Carlo sweep corresponds to Metropolis update of all 
I? nodes. We used a lagged Fibonacci pseudo-random number generator. For the 
simulations we used both scalar and parallel machines: a MASPAR MP1 massively 
parallel processor, a 12-node IBM SP2 and an 8-node DEC Alpha farm. The SP2 
and the Alphas were used as single independent CPUs, but we used a parallel Monte 
Carlo algorithm for the simulation of the largest lattice (L = 128) on the MP1. We 
show in Table p] the amount of thermalised data collected for the various lattice sizes 
and couplings in the flat phase. For the largest lattice we thermalised our surfaces 
by discarding about 10 7 sweeps. In addition we have performed simulations close to 




(5) 






(7) 
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L 


K = 1.1 


k = 2.0 


TR ~ L* 


32 


31 x 10 s 


26 x 10 s 


3 x 10 4 


46 


51 x 10 6 


42 x 10 6 


7 x 10 4 


64 


47 x 10 6 


44 x 10 6 


1.2 x 10 5 


128 


74 x 10 6 




5 x 10 5 



Table 1: The number of thermalised sweeps collected per data point in the flat 
phase. The last column indicates the autocorrelation time for the slowest mode 
in the system, the radius of gyration. The autocorrelation time is comparable for 
both values of k. The critical slowing down exponent z ~ 2, as expected for a local 
algorithm. 



the crumpling transition. This work is still in progress but preliminary results are 
discussed in Section [j| 

Previous studies of the critical behaviour of crystalline surfaces used more elab- 
orate simulation algorithms which combine a Langevin update with Fast Fourier 
Acceleration [15-^]. This algorithm is known to be more effective in reducing criti- 
cal slowing down. We chose not to use this approach for several reasons. First of all 
it is very hard to implement the Fast Fourier Transform (FFT) on a 2-dimensional 
surface with free boundaries. Secondly the parallel computer we are using is less 
efficient on floating point intensive algorithms, such as Fourier Transforms. Finally, 
the Langevin algorithm suffers from systematic errors induced by the finite time 
step At. This necessitates an extrapolation to At = 0, which can itself become time 
consuming. 

Our statistical errors were estimated using two methods: the first is a direct 
measurement of the autocorrelation in the data and a corresponding correction of the 
standard deviation. The second is the standard jackknife algorithm. Both methods 
give consistent results. 



2.2 Boundary Effects 

The optimal shape for a triangulated surface with free boundaries is a hexagon. In 
our simulations the use of the parallel computer MP1 makes the hexagonal mesh 
inconvenient, since the layout of the CPUs is a square grid. When we map the 
regular triangulation to the square grid, the resulting surface has a rhomboidal 
shape, as can be seen from Fig. |2[ In particular the regions in the shaded areas of 
Fig. § will be strongly influenced by the boundaries. For a generic observable O a 
we want to be able to quantify the effect of the boundary and of the anisotropy. In 
order to achieve this we integrate the observable over hexagonal subsets of the mesh, 
centred with respect to the surface. Fig. |2| shows two of these subsets with darker 
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Figure 2: Two of the hexagons used to define the averages of the integrated observ- 
ables. The arrow indicates the radius of the largest hexagon. 



lines. For a surface of linear size L we construct L/2 such integrated observables by 

°r = E °° (8) 

where H r is a hexagon of radius r and N r the appropriate normalisation. The shaded 
areas of Fig. ^ are discarded from the integration. By looking at larger and larger 
hexagons we can see when the boundaries start to affect the data. For very small 
hexagons the discretisation effects are large. In practice we find that the results are 
strongly influenced by the boundary for hexagons of radius r > L/4. 

For non-integrated observables, such as the normal-normal correlation function, 
translational invariance in the surface is broken by the presence of the free boundary. 
Thus we always consider the correlation function from the centre of the surface to 
all the other nodes. The effect of the boundary can then be seen clearly on the 
correlations at a distance of order L/4. This boundary data is discarded from the 
fits. 



3 Crumpling Transition 

Before examining the flat phase of the model of Eq. (|j) we would like to review 
the existing evidence for a "crumpling" transition. In recent years the crumpling 
transition has been the focus of extensive numerical and analytical investigation. 
Within the condensed-matter community it is customary to work with models with 
effective potentials like Eq. (0) and free boundaries [10,22,23]. Numerical simulations 



of these models provide direct evidence for a phase transition, such as a diverging 
specific heat, although the accuracy is not yet sufficient for a reliable determination 



8 



0.2 0.3 0.4 0.5 0.6 J^- 

(1 + K) 

Figure 3: The specific heat versus bending rigidity. 



of the specific heat exponent a. When strong self-avoidance is included in the models 
described by Eq. (0), there is numerical evidence that the crumpling transition 
disappears altogether, with the system being flat for all bending rigidities 28]. 



Some studies of flexible impenetrable plaquettes, however, do find a crumpled phase 
- pHR . The Gaussian spring models have been studied numerically in fTl, 15-21, 32| 



using periodic boundary conditions, with emphasis on the precise nature of the phase 
transition. A growing peak in the specific heat is observed and the best estimate of 
the related critical exponent is a = 0.58(10) |2lj| . 

There is thus strong evidence that models of phantom polymerised membranes 
have a continuous phase transition. We are also currently investigating this transi- 
tion, and in the rest of this section we discuss our preliminary results. 

3.1 Specific Heat 

Let us now turn to the energy fluctuations in our model. We write the total bending 
energy as 

S e = ]T n a • n b . (9) 

(ab) 
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Denoting the number of links (or edges) in the surface by N e , it is simple to show JT 
that the specific heat is given by 

Cv = 3 -^ Y ^ + ^({S!)-(S e f). (10) 

Here the brackets indicate a statistical average over surfaces^]. We henceforth drop 
the constant piece from our analysis. 

We plot in Fig. || the measured specific heat versus bending rigidity, for surfaces 
consisting of up to 64 2 nodes. As expected we see a sharp peak at tz ~ 0.79 growing 
with system size. The critical behaviour of Cy close to the phase transition is 
governed by an exponent a, Cy ~ \k — K c \~ a , and for a < 1 the phase transition 
is continuous (as the first derivative of the free energy does not diverge). Hence it 
is important to determine the value of a. The most convenient way of doing so is 
using finite size scaling, which predicts that the value of the peak should scale with 
volume as 

Cy « C + Cl L", (11) 

where, assuming hyperscaling, the specific heat exponent a = and cq and c\ 
are non-universal constants. Our best estimate of u>, from the data shown in Fig. || 
is oj = 0.5(1), consistent with previous results [21, [3^]. The corresponding value of 
a is 0.4(1). 

3.2 Radius of Gyration 

While the specific heat peak is one signal for the existence of a phase transition, 
it reveals little about the nature of the phases each side of the transition. An 
observable more sensitive to the geometry of the surface is the square of the radius 
of gyration, 

which measures the average spatial extent of the surface. R g defines a linear length 
scale for the surface and can be used to define a fractal or Hausdorff dimension dn 
in the embedding space via 

R g ~N 1/dH . (13) 

The Hausdorff dimension dn is related to the conventional Flory exponent v (R g ~ 
L u ) via dn = 2/v. For a flat surface, k > n c , the radius of gyration scales linearly 
with the internal size, and hence dn = 2. In the high-temperature phase, on the 
other hand, R g scales logarithmically with the volume of the surface, R g ~ log(iV). 
In this case we say that the Hausdorff dimension is infinite (y = 0). This justifies the 
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As our statistics close to the phase transition are not as good as in the flat phase, we have not 



attempted to use the more sophisticated method described in Section 2.2 to eliminate boundary 
effects from this data. 
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K 

Figure 4: The radius of gyration as a function of k for various system sizes. In the 
inset we show a power-law fit for R g vs. L at k = 1.1 for L up to 128. We get 
v = 0.95(5). 



terminology crumpled phase. This behavior can be computed exactly for k = [11| 
while mean field theory or numerical methods are necessary for k < k c [|l5|-21 
32 1 . Experimentally one can determine the Hausdorff dimension by measuring the 
structure function of diffracted light. A comparison with numerical simulations can 
be found in Ref. [33]. 

At the transition itself one might expect an intermediate behaviour (semi- 
crumpled surfaces) with dn > 2. Indeed this has been observed in [pj}| , where 
it is claimed that dn = 4 at n c . 

In Fig. |] we show our measurements of the radius of gyration versus the bending 
rigidity for surfaces up to 64 2 nodes. As expected, we see a dramatic change in their 
spatial extent as we pass through the phase transition. The surfaces literally blow 
up. This is better illustrated in Fig. || where we show snapshots of the surfaces: a 
in the crumpled phase, b at the phase transition and c in the flat phase. 

We have only determined directly the Hausdorff dimension in the flat phase 
(for k = 1.1), where we have good statistics for surfaces up to 128 2 nodes. The 
corresponding scaling plot is included in Fig. ||. A linear fit for surfaces larger than 
16 2 yields d H = 2.1(1), or v = 0.95(5), as expected for a flat surface. 
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Figure 5: Snapshot of the surfaces with L = 46 at (a) k = 0.5, (b) k = 0.8 and (c) 
k = 2.0. The average link length in each surface is comparable. 



3.3 Normal-Normal Correlation Functions 

As mentioned in Section |], the different phases of a crystalline membrane can also 
be distinguished by the behaviour of the surface normals. In the flat phase the 
normals will have true long-range order, with the correlation function approaching 
a non-zero asymptote. In the crumpled phase the normals eventually decorrelate. 

We define the normal-normal correlation function G(r) as the scalar product 
of two normals on the surface separated by a distance of r along the surface: 

G(r) = (n -n r ). (14) 

Here o refers to the centre of the surface. In the crumpled phase we expect G(r) 
to decay rapidly to zero. In fact in the exactly solvable k = Gaussian model one 
finds that 

where ci and C2 are cutoff-dependent constants and 6(r) is the two-dimensional 



regularised 5- function [11]. In the discrete case, the normals become decorrelated 
over a few lattice steps. 
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We have measured G(r) for several bending rigidities, both in the crumpled 
and flat phase. In Fig. ^ we illustrate this on lattices with 64 2 nodes. For a more 
detailed discussion of how we perform the measurements we refer to Appendix B. 

In the crumpled phase G(r) decays very rapidly to zero and, for k = 0, we indeed 
see a negative correlation between the normals at short distances, as expected from 
Eq. (|l5|). On a highly crumpled surface neighbouring triangles fold on each other. 
As k increases the normals become positively correlated at short distance, although 
G(r) still becomes negative (for k < k c ) before decaying to zero. 

At the critical point the normal-normal correlation function is expected to decay 
algebraically, 

G(r) ~ ±, (16) 

with an exponent fj different from the Gaussian value 4. This is consistent with 
our numerical results. A simple scaling argument ]23[ shows that fj is related to the 
size (Flory) exponent v by fj = 4(1 — v). As we enter the critical region (k ~ 0.79) 
the normal-normal correlation function still decays to zero, but now stays positive 
for all values of r. It is also clear that it decays to zero more slowly than in the 
crumpled phase. A crude estimate of the exponent fj, using a surface with 64 2 nodes 
and k = 0.79, yields fj ~ 0.71(5). This implies a size exponent v = 0.82(1) and 
hence Hausdorff dimension du = 2.4(2). This is consistent with previous numerical 
results |^3| and maybe be compared with two theoretical predictions: du = 3 from 
a 1/d renormalisation group calculation applied to d = 3, and du = 2.73 from 
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the self-consistent screening approximation [35|. 

Finally, we show in Fig. || a measurement of G(r) for n = 1.1, i.e. in the 
flat phase. In that case we see a non-zero asymptote indicating true long-range- 
order in the normals. Fitting G(r) to an algebraic decay plus a constant term 
excludes convincingly the possibility of a slow fall-off to zero, although eventually 
G(r) becomes zero due to boundary effects. This behaviour is found consistently 
for surfaces of various sizes. We will return to the behaviour of the normal-normal 
correlation function in the flat phase in Section 4.4. 



4 Flat Phase 

Here we shall describe the current theoretical understanding of the flat phase of 
crystalline membranes. Given the Mermin- Wagner theorem it is important to un- 
derstand what could stabilise an ordered phase in this two-dimensional system. This 
is most easily described in the Monge representation of a surface 

x CT = h a z + r a , (17) 

where h a is the height of the surface w.r.t. the base plane x — y and r a is the 
projection of x a on the base plane. We define the phonon fluctuations u ff of the 
surface by 

IV = s CT + u CT (18) 

where s CT are the equilibrium positions of the nodes. An effective Hamiltonian for 
the flat phase [§,|3(J is a sum of bending and elastic stretching energies 

H{h, u) = | J d 2 a (V 2 hf + 1 J d 2 a{2n u 2 a/3 + A u^), (19) 

where /i and A are the bare in-plane Lame coefficients, k is the bending rigidity and 
u a p is the strain tensor. This tensor measures the deformation of the induced metric 
from the flat metric, and is given by 

u a p = ^ (d« x • d/?x - <W0 (20) 
= \{^ a up + Vpu a + V a hV p h}. (21) 

to linear order in u. The indices a, run over the internal coordinate a = (a%, 02). 

One can integrate out the phonon degrees of freedom by performing the Gaus- 
sian integration |8|. One finds that the phonons give rise to an effective long-range 
two-point interaction for the scalar curvature. This interaction flattens the surface 
by stiffening the bending rigidity at large distances. 

A more precise understanding may be obtained from self-consistent calcula- 
tions of the renormalised bending rigidity || ^] , mean- field calculations |37j , the 
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e-expansion (e = 4 — D, with D the dimensionality of the surface) [38, 3^] and 
the 1/d expansion (where d is the embedding dimension) j34|,|I|. Generally one 
finds that the renormalised bending rigidity is scale dependent with a non- vanishing 
anomalous dimension 

K R [q) ~ q~\ (22) 

At long wavelength Kji(q) is stiffened by short wavelength undulations of the mem- 
brane. It is also found that the elastic moduli [i and A are softened at long wavelength 



(23) 



The exponents r\ and r\ u are not independent. Rotational invariance [38], or self 
consistent integral equations for the renormalised bending rigidity, imply 



V 



1 _Vu 
2 ' 



(24) 



An immediate implication of these anomalous dimensions or power-law singularities 
are non-trivial roughness exponents and height and in-plane displacement fluctua- 
tions or correlations. Defining a roughness exponent £ by the finite-size-scaling of 
the mean-square height fluctuations in a box of size L 



(h 2 



we see from Eq. (|l^) that 



(h 2 



d 2 q 



(2TT) 2 K R (q)q 4 



d 2 q 1 
(2tt) 2 q^ 



L 2 -\ 



(25) 



(26) 



implying £ = 1 — i]/2. In the above a is a short-distance regularisation cut-off. 
Similarly there is a non-trivial exponent for phonon fluctuations [41| 



(|u 2 |) 



L r,u 



(27) 



In the framework of the self-consistent screening approximation it is possible to ob- 
tain analytic predictions for these exponents [35]. By assuming the scaling relations 
of Eqs. (22) and (23) one can sum the terms in the perturbative expansion which 
renormalise Kn(q) and solve for the exponent rj. In subsequent sections we shall 
describe our numerical measurements of these exponents. 

Finally, a measure of long-range order in the flat phase is provided by the 
normal-normal correlation function. In the Monge representation a normal to the 
surface at point a is given by: 



(-d ai h,-d a2 h, 1) 



(28) 
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Figure 7: The distribution of eigenvalues p(X) of the inertia tensor. The number of 
nodes is 46 2 and the hexagon radius is 17. 

where a\ and 02 are the components of the intrinsic coordinate a. To compute 
(n CT • n D ) we rotate the surface so that the normal at the origin coincides with the z 
axis. Hence 



n CT • (0, 0, 1) 



1 



1 



1 



IV/U 



The correlation function then follows from the height-field propagator, 
( n<T . no ) ~ 1 - - / ^ -^-j « C+ - , 



2 7l (2vr) 2 /c fl (g) 



(29) 



(30) 



where r is the geodesic distance between the point at a and the origin o in the 
embedding space. Thus a power-law decay of the normal-normal correlation function 
to a non-zero asymptote provides a direct measure of the bending rigidity anomalous 
exponent 77. 



4.1 Shape Tensor 

We now return to our numerical results. Detailed information about the shape 
of the surface in the embedding space is provided by examining averages of the 
eigenvalues of the inertia tensor. More precisely we study the shape tensor, defined 
as the anisotropic part of the inertia tensor: 

1 

3iV 



(31) 
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where i and j refer to the components of x. The full inertia tensor has an additional 
isotropic contribution proportional to dij — this scales like Rg. As the shape tensor 
Iij is strongly influenced by the boundary (in fact contributions from the boundary 
region dominate the sum Eq. (31)), we restrict our measurements to a hexagonal 
subset of the mesh, as described in Section 2.2 . If we refer the surface to the body- 
fixed frame defined by the eigenvectors and its centre of mass, the eigenvalues of 
Iij, Xi, are nothing but the dispersion of the i-th component of x, averaged over the 
surface 

*i = (-5*T E (32) 




where N r is the number of nodes inside a hexagon of radius r. 

We obtain the eigenvalues Aj by diagonalising ly. The distribution of these 
eigenvalues, p(A), is distinctly different in the two phases. In Fig. ^ we show three 
examples of this. In the crumpled phase (a) p(X) has a single peak, implying a three- 
fold degeneracy of the eigenvalues^. This is a simple consequence of the isotropy 
of a crumpled surface in the embedding space. The position of the peak indicates 
the average extension of the surface, while its width is related to the fluctuations 
about that value. At the phase transition (6) there is still a single peak in p(X), but 
now accompanied by a significant tail extending to large eigenvalues. This is due to 
increasing fluctuations in the size of the surface in the critical region. 

The behaviour is very different for a flat surface (c). We now see two well 
resolved peaks in p(X) indicating, as expected, a broken 0(3) symmetry. There is a 
single minimum eigenvalue, corresponding to the left peak, which can be identified 
with the average square thickness of the surface. But there are also two almost 
degenerate large eigenvalues, as can be established by measuring the area of the 
right hand peak. This is a result of the remnant 0(2) symmetry in the plane. If 
we did not restrict our measurements to a hexagonal subset of the mesh, we would 
actually see three peaks in p(A), due to the anisotropy of the boundary. 



4.2 Roughness Exponent 

To measure the height fluctuations in the flat phase, and the corresponding rough- 
ness exponent (, we need an estimate of the out-of-plane fluctuations. This is 
provided by the minimum eigenvalue of the shape tensor which, as discussed in last 
section, is simply the average squared height. Hence 

Xmin ~ Q d 2 ahlj ~ L 2 «. (33) 

In Fig. [8| we plot the minimum eigenvalue A m j n versus the (normalised) radius of 
the hexagonal subsets we used, rjj = 2r/L. This is shown for four different size 

4 In fact this degeneracy is not likely to be exact. In a body-fixed frame there is always a hierarchy 
of eigenvalues. Fig. ma shows that this effect is very small. 
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Figure 8: The minimum eigenvalue X m in versus the radius of the hexagon for various 
lattice sizes and k = 1.1. The solid lines are obtained from a polynomial fit and the 
points are simply to guide the eye. 



surfaces at k = 1.1. A comment on how we treat the data: as we can only use 
hexagonal subsets of certain radii (discrete units of the lattice spacing), measure- 
ments on surfaces of different sizes will yield measurements at different values of the 
normalised radius r#. To compare measurements from surfaces of different sizes, at 
the same value of r#, we must therefore interpolate between the data points. These 
are the solid lines in Fig. ||. For the interpolation we used a polynomial fit, with the 
degree of the polynomials sufficiently large to yield a stable fit. We checked that 
different interpolation methods did not affect the results. The roughness exponent £ 
is then determined from the finite-size scaling of the minimum eigenvalue for a fixed 
value of ru. The result is shown in Fig. ^. The solid line is C( r H) while the dashed 
lines indicate the error. There are large discretisation effects for small values of r# , 
as expected. For 0.2 < th < 0.5 there is a reasonably stable value of the rough- 
ness exponent. Our best estimate from this intermediate region is C, = 0.64(2). This 
should be compared to the theoretical predictions C = -590 pH] and £ = 2/3 [|4[] and 
from measurements on the spectrin network £ = 0.65(10) [||. Previous numerical 
investigations of models of tethered surfaces have found a wide range of values for 
C- These include 0.5 fH, 0.53 @, 0.56(2) f|, 0.6 

. From the scaling relation 77 = 2(1 



0.64 m,|4|,|4|, °' 65 
C) we find r] = 0.72(4). For 

larger values of th we see clear evidence of boundary effects. Indeed, if we did not 



and 0.70 31 
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Figure 9: The roughness exponent £ versus the normalised radius of the hexagon. 
We extract the value of £ from the plateau region. The dashed lines indicate the 
size of the error bars. The bending rigidity is k = 1.1. 



work with hexagonal subsets, we would not be able to extract a reliable estimate 
for (. 

We have also measured C( r H) for k = 2.0, although the largest surface simulated 
in this case had only 64 2 nodes. Once again a finite-size-scaling analysis yields a 
stable value of £ in the same interval of r#. In this case the result is C, = 0.71(2). 
This value is 3cr larger than the corresponding value at the same lattice size for 
k = 1.1. We believe that this is due to larger finite size effects and that our results 
are consistent with having universal critical exponents in the flat phase as predicted 
by the fixed point of Refs. ||,||,|9|,|^. 



4.3 Phonon Fluctuations 

We now examine the issue of the in-plane elastic degrees of freedom in the flat 
phase of our model. For convenience, we rotate the surface so that the eigenvector 
associated to the smallest eigenvalue points in the z direction. Projecting the surface 
onto the x—y plane gives us a discretised analogue of the field r CT of Eq. (|l7|). 

The first step in the analysis must be to determine the phonon field u ff . We 
therefore need to determine the field s CT giving the equilibrium positions of the nodes. 



As before, we restrict our analysis to the hexagonal subsets H r of Section 12^. We 
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Figure 10: The exponent function of the hexagon radius. The dashed lines 

indicate the size of the error bars. 



assume that the equilibrium positions of the nodes lie exactly on a regular hexagon 
in the x— y plane. In the course of the Monte Carlo simulation, the surface fluctuates 
in the embedding space so that the orientation of its principal axes and its overall 
extent change constantly. Thus we need to find the regular hexagon s CT for each 
configuration we analyse. We define the functional 

F = Edv-s.) 2 = £ u 2, (34) 

a a 

and we choose the equilibrium position of the mesh to be best represented by the 
hexagon which minimises F. The regular hexagon s can be parametrised by the 
position of its centre, a rotation angle 9 and a scale factor £. The centre of the 
hexagon is trivially set to coincide with the centre of mass of the projected surface. 
Minimizing the functional F with respect to the angle 9 and the scale £ yields a 
unique solution up to the six- fold symmetry of the regular hexagon 9 — ► 9 + ir/6. 
This six-fold degeneracy of the minimum of F can be eliminated by requiring that 
the internal labels of the s a overlap with the ones of the r CT . 



We determine the characteristic scaling exponent rj u , defined in Eq. |27|, by a 



power law fit of (u ) rH to L, for fixed values of ru- Fig. 10 shows the result of the 



fit for k = 1.1. We extract the value of the exponent from the plateau in the figure. 
This gives rj u = 0.50(1) at k = 1.1. The corresponding value of r\ obtained from the 
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G(r) 




Figure 11: The G(r) fall-off for L = 32, 46, 64 and 128. Representative data points 
are shown to guide the eye. The dashed lines are the best fits to Eq. (pOl). 



scaling relation of Eq. fl24|) is rj = 0.750(5). We also quote the results of the same 
analysis at K = 2.0: rj u = 0.40(1). 



4.4 Correlations in the Flat Phase 

In this section we treat the normal-normal correlation function in the flat phase. 
The measured correlation functions are shown in Fi g. p4| . We fit this correlation 
function to the expected power-law behaviour of Eq. (|30|) 



G(r) = C + 



B 



(35) 



using a correlated least squares algorithm. We find a non-zero asymptote C whose 
value tends to increase with lattice size. This is shown in Fig. |l^. The correlated 
least squares algorithm attempts to minimise a chi-squared function of the form 



N 



X 



N(N - P) 



J2(yi - y*{Pi)) c ij(yj - y*j{pi))- 



(36) 



The data is fitted to a functional form y*(pi) where i is the (lattice) geodesic distance 
and the quantities {pi},i = 1 . . . P, are the variational parameters in the fit. The 



matrix Cjj is the inverse of the correlation matrix 



(M-i) 



IJ 



(37) 
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0.2 ~i ' 1 1 ' ' 1 ' ' 1 1 ' 1 1 ' 1 1 ' ' 1 r 

0.2 0.4 0.6 0.8 

r/L 

Figure 12: A plot of G(r) with r rescaled to show the finite-size effects on the 
asymptote (k = 1.1). This rescaling demonstrates that the asymptote stabilizes as 
the lattice size increases. 



where 

M rr , = (G{r)G{r')) c . (38) 

Note that for uncorrelated data the matrix M is diagonal, and Eq. (^6|) reduces 
to the usual definiition of x 2 . The inversion of this matrix is typically a delicate 
operation. Because of limited statistics it will often be close to singular and the 
smallest eigenvalues will only be poorly estimated from the data. In this situation 
it is necessary, and indeed correct, to restrict the inversion to an appropriate sub- 
space that is spanned by eigenvectors with eigenvalues which are large enough to be 
estimated reliably from the data. This can be achieved through singular value de- 
composition techniques. The dimension of the subspace is referred to as the singular 
value decomposition (SVD) cut. For more details we refer the reader to [51]. 

In assessing the results of the fitting procedure we examined cases with a range 
of values of initial and final distances and SVD cuts. To obtain a x 2 °f order unity 
we usually had to discard data from the last third of the path to the boundary and 
roughly one quarter of the eigenvalues. 

In Table || we show the fitted values of the asymptote C and the exponent r\. 
The fits are obtained including all the short distance data and we show only errors 
from the fit. It is clear that a more important source of error is finite size effects. 
Clearly both the constant C and the exponent r\ exhibit a shift with L well in excess 
of the fitted errors. But this shift is systematic and indeed can be quantified by 
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L 


K = 1.1 

C j) x 2 /dof 


k = 2.0 
C f] X 2 /dof 


32 
46 
64 
128 


0.181(16) 0.331(7) 8.6 
0.274(9) 0.397(5) 1.66 
0.321(5) 0.447(6) 0.96 
0.383(6) 0.521(6) 1.14 


0.240(32) 0.141(5) 14.8 
0.309(21) 0.154(4) 2.32 
0.448(30) 0.203(4) 2.19 



Table 2: Fit to the correlation function G{r). 



Train 


The constant C 
L = 46 L = 64 L = 128 


The exponent rj 
L = 46 L = 64 L = 128 


1 

2 
3 


0.274(9) 0.321(5) 0.383(6) 
0.262(15) 0.325(7) 0.388(7) 
0.262(17) 0.327(7) 0.405(7) 


0.397(5) 0.458(7) 0.521(6) 
0.376(10) 0.492(9) 0.546(9) 
0.388(14) 0.481(9) 0.658(15) 



Table 3: Effect of the short distance data (k = 1.1). 



assuming a naive 1/L correction term to the correlation function G(r). In that case 
the leading correction to C will be 1/L, while it is log(L)/L for the exponent r\. 
This conjecture is in excellent agreement with our data, at least for k = 1.1, and 
implies infinite volume values of C ~ 0.45 and r] ~ 0.62. This is in quantitative 
agreement with our previous estimates. 

Another source of systematic error stems from discarding some of the short 
distance data. In Table |3| we show the results of fits to C and r] obtained by 
discarding data at distance r < r m i n . It is clear that the values change with r m i n , 
although not a great deal. Discarding short distance data increases the value of r] 
for the large lattices. This improves the comparison with the expected value. 

Finally, a few comments on the data at k = 2.0. The values of r] quoted in 
Table || do not agree well with the values obtained at k = 1.1. But it should be 
mentioned that the fits to Eq. (35) are not as good in this case, which is reflected 
in higher values of the \ 2 ■ 



5 Conclusions 

In this paper we used a large scale Monte Carlo simulation to show that an extremely 
simple model of a crystalline membrane has a well-defined flat phase. In this phase 
the critical exponents are consistent both with previous simulations of tethered 
membranes and with theoretical predictions. In particular, the model we study has 
dynamically generated elastic moduli. The flat phase is convincingly demonstrated 
by the existence of a non-zero asymptote for the normal-normal correlation function, 
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which strengthens with increasing system size. The flat phase exponents we find at 
k = 1.1 are: size (Flory) exponent v = 0.95(5) (dn = 2.1(1)), roughness exponent 
C = 0.64(2), rj u = 0.50(1). A check on the consistency of our results for £ and rj u 
is obtained by comparing their respective scaling predictions for r\. Our value of £ 
implies 77 ~ 0.72 and our value of rj u implies rj ~ 0.75. These are consistent within 
our statistical errors. A direct measurement of r\ from the power law decay of the 
normal-normal correlation functions is more difficult and is discussed in Section <L4 . 
For L = 128 and 

r min — 3 we obtain rj — 0.658(15). At higher bending rigidity we 
obtain somewhat different exponents, but we believe this is due to larger finite-size 
effects. 

For completeness, we also establish that the model we consider has a crumpling 
transition. Preliminary determination of the critical exponents at the transition 
gives results consistent with existing simulations for related models. 
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Appendix A: Elastic Constants in the Gaussian Model 

As we noted in the introduction, there are several potential pathologies in a model 
in which the tethering potential contains no hard core repulsion, such as the one we 
treat in this paper. The effective elastic constants may vanish or be too weak to 
generate a stable flat phase. Even if the model possesses a flat phase it may belong 
to a different universality class from a model with bare elastic constants. 

The simplest argument supporting such concerns arises from a mean field theory 
calculation of the elastic moduli of our discrete tethering potential along the lines 
of |33]]. Consider an equilibrium spring length parameter a in the pair potential of 
Eq. (|): 

H T [a) = ]T (|x a -x CT ,|-a) 2 (39) 
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^EE(i x -- x -'i-«) 2 - (40) 



a a' 



In the second sum a' runs over the neighbour nodes of a. Note that we have chosen 
units such that the spring constant k/^ksT) = 1 and x, u and a are dimensionless. 
For a sufficiently large we may describe the location of the nodes x CT by 

(x CT -x a ,) a = {8 aP + u a p) sj a , , (41) 

where s a spans a regular hexagonal lattice, s acr i = s a — s a i and \s a(T '\ = a. a and 
(3 = 1, 2 are intrinsic coordinate indices on the surface. Then 



X,t X^/ 



(8 aj3 + 2u a p + UayUfr) s^s^, (42) 



= + (2u Q/3 + Ua-yU/3-y) S^'S^, . (43) 

At this point we redefine u' Q( g = u a p + \u a ^u^ and subsequently drop the primes. 
This will not affect the end result to quadratic order. Expanding for small fluctua- 
tions Uap < 1 we have 

|x CT -x CT ,| = (a 2 + 2^8^,8^,) 1 (44) 

~ a + - u apS%a' s L> - ^S^a^S^S^S^S^, + . . . . (45) 



To evaluate Eq. (|40|), consider the 6 unit lattice vectors d;, of a regular hexagonal 
lattice where b = 1, . . . , 6. Performing the sum over the neighbours a' gives 



n T {a) = \Y,il( au ^ d % d b) 2 ( 4 6) 

<r 6=1 

= I q2 E (V 7 ^ 5 + rt^ 7 + rt 7 " 5 ) Ua(3 u l5 (47) 
= -a 2 (Zuapupa + U^j . (48) 

(T 

For sufficiently large N, and a fixed, one can approximate the above discrete sum 
with an integral. 



<H f \ ~ f ^ ^ 2 , \/3 2 2 \ 

H T {a)K, J — I —a u a pup a + —a « 77 I 



(49) 



where e is the intrinsic lattice spacing. Often one is interested in the case e ~ a. We 
keep the e and a dependence since we are interested in the case a = 0, e ^ 0. 
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To determine the elastic constants recall that the most general continuum action 
for a crystalline membrane is [37] 



n 



d 2 i 



7j<9 2 x • 5 2 x + u (<9 Q x • <9/?x) 2 + 



+ v (<9 7 x • d-yx) 2 + - (<9 a x • dpx) 



^<9 2 x • <9 2 x + fiu a/3 ui3 a + -it 2 7 + r« 77 



(50) 
(51) 



where u a p = \{d a x ■ c^x — 5 a p) is the strain tensor Eq. (p0|). The couplings in the 
two expressions in Eq. (|50|) are related by 



4 ' 



t = T — (JL — A. 



(52) 



Thus the bare elastic constants, dimensionless in our units, arising from the lattice 
action (^0|) are 



fj, = A 



\/3 a 



2 

4 e 2 " 



r = 0. 



(53) 



Notice that if e = a the elastic constants are independent of the lattice spacing. The 
model we have studied corresponds to the limit a — * and e fixed. In this limit 
fluctuations become large and the above calculation breaks down. 



Appendix B: Measuring Correlation Functions 

A necessary premise for measuring the normal-normal correlations function G(r) 
is that we know the distance between two triangles on the surface. This distance 
can either be measured in the intrinsic metric (in which case it is trivial) or in the 
induced metric. For comparison we tested both methods. 

Let us first describe our algorithm for determining distances in the induced 
metric at the discrete level. Given two points on the surface, we must find the 
shortest path, in the induced metric, along the surface. The problem here is that 
there are many definitions of "geodesic" which are equivalent in the continuum limit 
but differ at the discrete level. The algorithm we use is the following: given a triangle 
to we find the distance from its centre to the centre of all its neighbours. Then 
we find the distance from those triangles to their neighbours (excluding triangles 
already visited). Iterating this procedure we find the minimum distance to each 
triangle from to, subject to the constraint that we have to pass through the centre 
of each triangle traversed. This is a piecewise linear approximation to the shortest 
path, which should be good for sufficiently large paths. 

Similarly, for the intrinsic metric, we can either define the distance in units 
of jumps from triangle to triangle or, given that the vertices have explicit ((71,(72) 



26 



G(r) 



0.4 " 



0.7 " 



0.6 " 



0.5 " 




K= 1.1 



V =64 2 



0.3 







5 



10 



15 



20 



25 



30 



35 



40 



r 



Figure 13: The comparison of the three methods used to measure the geodesic 
distance r. Methods A and B use the intrinsic metric. For method A we construct 
a piecewise linear path through the centre of each triangle, while for method B 
we construct a straight line between the centre of the starting and ending triangle. 
Method C is the same as A, except that it uses the induced metric to compute the 
distances between the triangles. 

coordinates in the intrinsic formulation, we can define the shortest distance between 
them as a straight line. 

We have verified that these different definitions of distances give identical results 
for the normal-normal correlation function measured in the flat phase (modulo a 
trivial rescaling of the r-axis). This is illustrated in Fig. The relevance of this 
result goes beyond a simple consistency check: the fact that the geodesies defined 
intrinsically and extrinsically coincide means also that intrinsic and the extrinsic 
metric overlap in the flat phase. 

All of the normal-normal correlation functions presented in the paper are ob- 
tained using method A. 

Appendix C: Parallel Monte Carlo Algorithm 

In this appendix we describe the parallel algorithm used on the MASPAR MP1. This 
machine is an old-style massively parallel computer with 16384 CPUs arranged in a 
2 — D (128 x 128) mesh with nearest neighbour connections. Each individual CPU 
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Figure 14: The improved parallelisation scheme. On the left we show the connec- 
tivity of the 4 different surfaces in the simulation, while on the right we show the 
actual internal representation. It is clear that, on the right, no node ever has a 
neighbour of the same type. 

is a relatively small processor (8 bit) with no floating-point unit. 

A standard problem in using a parallel machine is the fact that the amount of 
parallelisation, and consequently the performance increase, is limited by the inter- 
dependence of the data. In order to ensure detailed balance in a Monte Carlo 
simulation, only a fraction of the lattice can be updated concurrently. In our case 
only 25% of the nodes can be updated. This translates into a huge performance loss 
as 75% of the nodes remain idle. 

In order to overcome this limitation we implemented an improved parallelisation 
scheme. Instead of simulating one surface we consider 4 independent Monte Carlo 
simulations. The 4 corresponding meshes are "interleaved" in 4 arrays which store 
the node positions. Each surface is distributed onto the 4 arrays as shown in Fig. 



14j. The parallel machine updates one array at a time, but now each array holds 
independent data and therefore all of it can be updated in parallel. 

We have compared the performance of this algorithm to the traditional paralleli- 
sation and to the sequential code. On the MASPAR, the traditional parallelisation, 
i.e. a Monte Carlo simulation of a single surface, can achieve a maximum speed of 
80 Mflops (millions of floating-point operations per second). Our improved code 
is capable of 280 Mflops, almost a four-fold increase. This number is to be com- 
pared with the peak performance of this machine which, measured with the Linpack 
method, is around 440 Mflops. Our scalar code on an IBM RS/6000 390 with a clock 
of 66.5 MHz has a speed of 17 Mflops, compared to a Linpack peak performance of 
54 Mflops. 
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The parallel algorithm also requires a careful implementation of the pseudo- 
random number generator. While the intrinsic MASPAR random number generator 
is known not to be reliable, using a sequential random number generator would be 
incredibly time consuming. The solution to this problem is to have an independent 
random number generator on each CPU, e.g. to regard it as an array-valued function. 
In order to avoid cross correlations between the random sequences generated by the 
parallel routine, a second standard random number routine is used to seed the array. 
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